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Abstract 

Consider two coupled oscillators, each modulating 
the other’s frequency. This system is governed by 
four parameters: the base frequency and modulation 
index for each oscillator. For some parameter values 
the system becomes unstable. The Lyapunov ex¬ 
ponent is used to measure the instability. Images of 
the parameter space are generated, with the number 
crunching implemented on graphics hardware using 
OpenGL. The mouse position over the displayed im¬ 
age is linked to realtime audio output, creating an 
audio-visual browser for the 4D parameter space. 
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1 Introduction 

Soft Rock EP [ClaudiusMaximus, 2005] and 
Soft Rock DVD [ClaudiusMaximus, 2006] ex¬ 
plored the transitions between order and chaos 
in coupled FM oscillators. A more recent con¬ 
tinuation of this project is to make a map of the 
parameter space of coupled FM oscillators on a 
perceptually relevant level and use it in live per¬ 
formance, choosing parameters on the basis of 
desired sound character. 

A bifurcation diagram produced by an ana¬ 
logue Moog synthesizer [Slater, 1998] and im¬ 
ages of Lyapunov fractals [Dewdney, 1991] were 
inspiration to apply the latter technique to the 
parameter space of coupled FM oscillators in 
the digital realm. 



Figure 1: Example output. 



2 Formulation 

2.1 Coupled FM Oscillators 
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Figure 2: Coupled FM oscillators in Pure-data. 

Consider the two coupled oscillators in Fig¬ 
ure 2. Pure-data’s model of interconnected 
components each with their own internal state 
maps poorly to GPU architecture. Consider¬ 
ing the whole system as one and flattening the 
internal state into a single phase space vector 
leads to the following formulation as a mutual 
recurrence relation: 

x n+1 = %(x n + I(/ x + m x cos(27ry n _ rf ))) 

Vn+l = %(Vn + I (fy + m y COs(27 TX n - d ))) 

where 

„ . . . . i . 440 t—69 

m = g^ 2 - 

Here x n , y n is the phase of each oscillator at 
time step n, d is a delay measured in sam¬ 
ples, f x , f y is the base frequency of each os¬ 
cillator as a MIDI note number, and m x , m y 
is the modulation index of each oscillator as a 
MIDI note number. %(t) performs wrapping 
into [0,1), with \t.\ being the flooring operation 
(the largest integer not greater than t). 

The four-dimensional parameter space vector 
will be written 

® = (fxi fyi Wlxi TTly) 

and the (2d + 2)-dimensional phase space vector 

2 = ( x n , y n , x n —i, y n ~ i, • • •, Xn—d, y n —d) 

with sample rate SR = 48000. For reasons ex¬ 
plained in Section 5.2, d = 1 will be fixed. 


2.2 Lyapunov Exponents 

Lyapunov exponents can be used to measure the 
stability (or otherwise) of a dynamical system. 
A good introduction is found in Chapter 4.3 
Lyapunov Exponent [Elert, 2007]. The defini¬ 
tion is covered in Chapter 13.7 Liapounov expo¬ 
nents and entropies [Falconer, 2003] which also 
relates it to measures of fractal dimension. 

The Lyapunov exponent A measures diver¬ 
gence in phase space: 
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An attracting orbit has A < 0 and a divergent 
(chaotic) orbit has A > 0. 

A modified norm is required to take into ac¬ 
count the wrapping of phase into [0,1): 
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For example the distance between 0.1 and 0.9 is 
properly 0.2 (not 0.8) because 0.1 can be phase- 
unwrapped to 1.1. 


2.3 Viewing Planes 


An image is 2D, which requires choosing a sub¬ 
set of the 4D parameter space to visualize. Two 
particular planes were chosen: 
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where (u, v ) is the coordinates of the pixel, ao 
is the centre of the view, and ro is the radius of 
the view. 

These planes were chosen because they are 
simple, while still being flexible enough to ex¬ 
plore the whole 4D space. The A + plane varies 
both oscillators in the same direction, while the 
A_ plane varies each oscillator in opposite di¬ 
rections. To center on a particular target point 
{fx,f y ,m x ,m y ) one might use the A + plane to 





center on the midpoint 

( fx + f y fx + f y m x + m y m x + m y \ 

V 2 ’ 2 ’ 2 ’ 2 ) 

and then switch to the A_ plane to break the 
(x, y) symmetry. 

3 Implementation 

The implementation uses OpenGL [Segal, 2013] 
and OpenGL Shading Language [Kessenich, 
2013] for computation and graphical rendering, 
GLUT [Kilgard, 1996] for windowing and input 
event handling, and JACK [Davis, 2013] for au¬ 
dio output. 

3.1 Introduction to Modern OpenGL 

Modern OpenGL has a programmable shader 
pipeline. Vertex attributes are read from vertex 
buffers and processed by vertex shaders. The 
outputs of the vertex shader (called varyings) 
are further manipulated by an optional geome¬ 
try shader stage. Geometry shaders can output 
a different vertex count to their input count, 
whereas vertex shaders are one-in one-out. The 
result of the geometry shader can be captured 
into another vertex buffer using transform feed¬ 
back. Following the geometry shader the prim¬ 
itives (points or triangles) are rasterized, and 
varyings interpolated across each primitive. Fi¬ 
nally a fragment shader takes these values and 
computes the colour at that pixel. The output 
of a fragment shader can be captured by attach¬ 
ing a texture to a framebuffer. 

3.2 Computation Overview 

To render an image a texture is first filled with 
(u, v ) coordinates using a framebuffer object 
and a fragment shader. This texture is copied to 
a vertex buffer object, interleaved with an initial 
phase space vector z = (0, 0, 0,0) and Lyapunov 
exponent statistics vector l = (0,0,0,0) for each 
point. 

Using a vertex shader, a is calculated from 
(u, v ) using Equation 3, and then a rough es¬ 
timate of the Lyapunov exponent is computed 
using Equation 2 by perturbing 2q(0) = Zo(0')+d 
with 5 small and performing t = 256 iterations 
of Equation 1. The first few repetitions are dis¬ 
carded, along with those resulting in —oo, and 
the rough A estimates are accumulated in l. 

Between each repetition the working set is 
compacted using a geometry shader. Points 
whose mean Lyapunov exponent estimate 
changed very little during the previous step are 


plotted and removed from the working set. The 
other points are kept to be refined further, di¬ 
recting the computational effort on the points 
that need it most: those slow to converge. 

To ensure user interface responsiveness, the 
computation is amortized over several frames. 
The target frame period is divided by the mea¬ 
sured time for one repetition to compute how 
many repetitions to perform that frame. The 
repetitions-per-frame increases as the working 
set becomes smaller. 

3.3 Noise Increases Stability 

At the end of each repetition z\ is kept instead 
of zq. This effectively adds a small amount 
of noise, counter-intuitively increasing stability. 
Noise allows more of the phase space to be ex¬ 
plored, and makes it more likely for the per¬ 
turbed orbit to reach an attracting part of the 
phase space. 

3.4 Dither Increases Quality 

To reduce grid sampling artifacts, (u, v ) is per¬ 
turbed within the bounds of its corresponding 
pixel before calculating the a parameter vector 
for each repetition. 

4 Results 

4.1 Examples 

Figure 3(a) shows the initial view on starting 
the interactive browser. Low frequencies to the 
left are stable even at high modulation index 
away from the central axis. High frequencies to 
the right become chaotic at progressively lower 
modulation index, (b) shows the A_ plane at 
the same location, (c) shows bands alternat¬ 
ing between stability and chaos. The bands be¬ 
come distorted and collapse as the modulation 
index and frequency increase, (d) shows its A_ 
plane, bands become rings. When the frequency 
is greatly increased, the shapes become more 
intricate, (e) exhibits spirals of stability, with 
similar spirals in the A_ plane in (f). 

When f x = f y and m x = rn y the A + plane has 
mirror symmetry about its horizontal axis, and 
the A _ plane has two-fold rotational symmetry 
about its centre. Breaking the symmetry and 
setting f x 7 ^ f y or m x / rn y leads to diverse 
forms. In particular Figure 3(h) has shapes that 
resemble those of Lyapunov space images of the 
logistic map. 

4.2 Interactive Explorer 

The implementation is an interactive audio¬ 
visual explorer for the parameter space of cou- 
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pled FM oscillators. Clicking with the mouse 
zooms the view about the clicked point. The left 
button (or scroll up) zooms in, the right button 
(or scroll down) zooms out, the middle button 
centers the view on the target point. Pressing 
the TAB key toggles between the A + and A_ 
planes in Equation 3, and Fll toggles full screen 
operation. 

While the GPU simulates and analyses one 
oscillator pair per pixel, the CPU simulates one 
oscillator pair with a determined from the pixel 
under the mouse pointer. The image acts as a 
map, a reference frame for chosing parameters 
to audition by moving the mouse. 

5 Conclusions 

5.1 Original Intent 

Earlier experiments used one Pure-data batch 
mode instance per CPU core each sending anal¬ 
ysis data to a realtime Pure-data instance. The 
analysis used various methods (including FFT 
for spectral statistics and the sigmund exter¬ 
nal for pitch tracking) to classify points into 
pitched (ordered, stable) or unpitched (chaotic, 
unstable) with measures of distortion or noisi¬ 
ness. Sadly this approach proved impractical as 
it achieved only tens of pixels per second, even 
with a fast multi-core CPU, and porting these 
signal analysis algorithms to massively-parallel 
programmable graphics hardware seemed to be 
too difficult. 

5.2 OpenGL Issues 

The current implementation is hardcoded with 
delay d = 1 and would be very awkward to 
generalize. OpenGL architecture limits each 
vertex attribute to four components with the 
maximum number of attributes typically lim¬ 
ited to sixteen. This totals 64 floats per ver¬ 
tex, 6 of which are needed for the pixel coor¬ 
dinates and Lyapunov exponent statistics accu¬ 
mulation. Therefore using OpenGL imposes a 
limit d < 28. For comparison the original ex¬ 
periments in Soft Rock EP used Pure-data’s de¬ 
fault block size of 64, with d = 32. Moreover, 
increasing d increases video memory consump¬ 
tion. With the maximum d = 27, browsing at 
1920 x 1080 resolution would require over 1GB. 

Future work on this project will look into 
using OpenCL, which provides a heterogenous 
CPU and GPU computation framework, in the 
hope that it will avoid the inherent awkward¬ 
ness of abusing OpenGL shaders to perform cal¬ 
culations. 


5.3 Audio Issues 

While the implementation works as intended, 
with d = 1 the sound is nowhere near as rich 
and varied as with d = 32. With small d there 
is much more very high frequency content in 
interesting-looking regions. There seem to be 
few if any regions of the parameter space with 
both interesting appearance and palatable au¬ 
dio frequencies at d = 1, while with high d there 
are parameters that generate sounds that fluc¬ 
tuate intermittently between smooth tones and 
noise. Visualization with high d has not been 
possible so far, so whether their neighbourhoods 
look as interesting as they sound remains an 
openquestion . 

Unfortunately, heavy use of the GPU in the 
interactive browser can block the operating sys¬ 
tem for too long and cause audible glitches 
(JACK xruns). This situation may change as 
free drivers continue to improve, allowing use of 
the browser in a live situation. 

5.4 Pretty Pictures 

Despite these shortcomings, I think the images 
look good. I plan to render a selection at high 
resolution and print postcards and posters. For 
huge images it is possible to divide the image 
plane into tiles and compute each tile in succes¬ 
sion, finally combining the pieces into one large 
picture. 

There is also scope for video work, moving 
and rotating the viewing plane through the 4D 
parameter space, with different shapes forming 
and collapsing over time. Rough benchmarks 
take 5-10 seconds per frame at 1920 x 1080, 
so it seems sensible to wait until faster cheaper 
graphics cards become available. 

6 Obtaining the Implementation 

The implementation was written on 
GNU/Linux Debian Wheezy running on a 
quad-core AMD64 processor with NVIDIA 
GTX 550Ti graphics card using propri¬ 
etary drivers. The source code is available at: 
https://gitorious.org/maximus/lyapunov-fm 
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